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Abstract 

The explicit stability constraint of the discontinu- 
ous Galerkin method applied to the diffusion operator 
decreases dramatically as the order of the method is 
increased. Block Jacobi and block Gauss-Seidel pre- 
conditioner operators are examined for their effective- 
ness at accelerating convergence. A Fourier analysis 
for methods of order 2 through 6 reveals that both 
preconditioner operators bound the eigenvalues of the 
discrete spatial operator. Additionally, in one dimen- 
sion, the eigenvalues are grouped into two or three 
regions that are invariant with order of the method. 
Local relaxation methods are constructed that rapidly 
damp high frequencies for arbitrarily large time step. 

Introduction 

The discontinuous Galerkin method allows compact 
spatial discretizations of any order to be formulated 
on arbitrary meshes. The compact discretizations are 
well suited to explicit time-marching methods for un- 
steady simulations, and are easily and efficiently im- 
plemented to run in a parallel environment. 1 * When 
applied to propagation and advect.ion equations, the 
stability constraint decreases only moderately - as the 
order of the method is increased. However, when ap- 
plied to the diffusion operator, the explicit stability 
limit decreases dramatically as the order of t he method 
is increased, 3 and some degree of implicitness is nec- 
essary. 

Implicit time marching methods ranging from back- 
ward Euler to a high-order implicit Runge-Kutta could 
be employed, depending on the level of temporal accu- 
racy required. Regardless of the approach, a globally 
implicit solution is required. The direct approach, a 
direct solver, destroys the compactness inherent to the 
discontinuous Galerkin method along with the associ- 
ated advantages. Thus we are motivated to focus on 

•Senior member, Senior Research Scientist, Computational 
Modeling and Simulation Branch 

t Professor, Division of Applied Mathematics 
Copyright © 2001 by the American Institute of Aeronautics 
and Astronautics, Inc. No copyright is asserted in the United 
States under Title 17, U.S. Code. The U.S. Government, has 
a royalty-free license to exercise all rights under the copyright 
claimed herein for Governmental Purposes. All other rights are 
reserved by the copyright owner. 


iterative methods that are local to the element and 
preserve the compact character of the discontinuous 
Galerkin spatial discretization. One obvious choice is 
a block Jacobi relaxation scheme in which some or all 
of the local contributions to the spatial operator are 
treated implicitly. 

This paper presents the Fourier analysis of block Ja- 
cobi and block Gauss-Seidel preconditioning operators 
for the discontinuous Galerkin method applied to the 
diffusion operator. Data from the analysis is used in 
the development of optimal relaxation methods that 
rapidly damp the high frequencies, and are well suited 
for use with multi-grid 7 techniques. Numerical tests 
are presented that verify the analysis. 

The first, section describes the discontinuous 
Galerkin method as applied to diffusion, introduces 
required notation and formulates the particular im- 
plicit problem that is the focus of this paper. I he 
second section presents a detailed Fourier analysis, in 
one dimension, of the preconditioning operators and 
the findings of t hat analysis. I he t hird section for- 
mulates several optimal relaxation methods based on 
data from the analysis, and demonstrates them in nu- 
merical tests. The fourth section extends the analysis 
to two dimensions and formulates optimal relaxation 
methods that are demonstrated in numerical tests. 

Methodology 

In this section, we describe in general terms the ap- 
plication of the discontinuous Galerkin method to a 
model diffusion equation of the form 

^ - V • [//Vu] = 0. (i) 

dt 

Although there are several ways in which the discontin- 
uous Galerkin method may be applied to a diffusion 
operator, 3 " 5 this analysis will focus on an approach 
that requires a first order equation. To this end, equa- 
tion (1) is recast, as a set of first, order equations: 


if-ftVu = 0. (2) 
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Spatial Discretization 

The discontinuous Galerkin method is applied on a 
discrete element, Qi obtained by partitioning the do- 
main into arbitrarily shaped elements. The solution 
in each element is approximated in terms of a set of 
local basis functions as u ft* m — 5Z{f} 

and q w qi = £ {l} 111 The g overnin 8 e( l ua ‘ 

tions are projected onto each member of the basis set 
to give a set of equations in each element that govern 
the evolution of the set of unknown variables {«/,,-}• 



dQ 


J [qi - Vui] dQ 
n 


0 

0. 


(3) 


The gradient term is eliminated by integrating by parts 
to obtain the weak form 


/[•«£ 

n L 

/ 


- V6i,i • qi 


dQ + ^ j'bi t ih(qii<fk) -ds = 0 

{fc} 


[6/^7} — V6/ t i • Ui] dQ -f 


e A- 


h(ii(,Uk)ds = 0 (4) 


where ft denotes a numerical flux, the subscript Ar de- 
notes the index of a neighboring element. dQi tk denotes 
the boundary between Qi and Q k , and ds is the out- 
ward boundary normal. Note that the variables {qi} 
are simply intermediate values that can be computed 
directly and explicitly from the solution. In the case of 
advection and propagation, the numerical flux is usu- 
ally an approximate Riemann flux. For diffusion, the 
flux is simply a weighted average of u or <f such that 


h(qi,qk) — (1 — 4>l,k )<f/ T 

h(ui , Uk) = + (1 - <f>i,k)uk, 

and (pi k - l “ Evaluation of the integrals and 

fluxes of equation (4) leads to a semi-discrete evolution 
equation of the form 


B 


= [A./ = - 



C = [<■; ,] = -M _ 1 I Jbt'ibkjds , 

Un,, fc 

and 

M = [m;j] = ^ bijbijdQ . 

To facilitate the analysis, it is assumed that the do- 
main has been uniformly partitioned such that the 
matrices A, B^ and are the same for all elements. 
Further more, in this study we will examine only the 
so called "mixed LR V scheme 3 ’ 5 in which <pi, k is either 
zero or one on each boundary segment and 

is assigned in a regular and repeating pattern. 

Temporal Discretization 

In previous works dealing with propagation,*"'" ex- 
plicit Runge-Kutta was used to advance the solution 
in time. However, the stability bound of such an 
explicit method is inversely proportional to the maxi- 
mum eigenvalue of the spatial operator. When discon- 
tinuous Galerkin is applied to the diffusion operator, 
the maximum eigenvalue grows in magnitude approx- 
imately proportional to p 4 , and the explicit approach 
becomes impractical for most problems. An implicit 
time marching method would clearly remedy the prob- 
lem; however, a direct inversion of the resulting mat rix 
would destroy the compactness of the discontinuous 
Galerkin. Thus in the present work, we consider im- 
plicit evolution methods that are solved by local relax- 
ation methods, and we are interested in the analysis 
of the relaxation process. 

Given a general evolution equation of the form 


du 

!h 


= fl(n), 


where R(u) denotes a spatial operator, all implicit 
methods involve solving an equation of the form 


= A . Q| + [(1 - ■ Q i + 4i,kC k • Q k] 

dt U-} 


Q, = AU, + + (1 - ^.*)C fc Ufc] , 

{*} 

where 



" Ul, l 




Ul o 

Q( = 

qi, 2 

u, = 

Ul, 3 

qi.3 


: „ 


_ 


A = [Sij] = M- 1 


f V6i,i ■ 6/ jdtt 


(5) 



( Ar ) 2 


R(u n - l ),R(u n - : *),...), 


(6) 


where the superscript n may denote either the time 
step (t = n A t) or a step within the Runge-Kutta 
sequence, and a and H depend on the specifics of the 
method. In the present analysis, we are concerned 
only in the construction and stability of an iterative 
method for the solution u n in equation (6). Thus, we 
assume the discrete evolution equation (6) is stable in 
time, possesses the desired order properties, and that 
old values of the solution u M l , u 7) , are known. 

Within the scope of a Fourier analysis, the terms 
involving the old values of the solution have no offec t 
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on the stability of a relaxation scheme. Thus it is suf- 
ficient to examine a much simplified form of equation 
(6) given by 


1l( A,v) = Aft(t>) - v = 0, (7) 

where A = a A t/( Ax) 2 , and v = u n . Including a in 
the A term eliminates all explicit dependence on the 
specific evolution method, and the conclusions of the 
following analysis apply to all time evolution meth- 
ods that can be written in the form of equation (6). 
Because this class of evolution methods includes first- 
order backward Euler, the following analysis also has 
implications for steady-state solution methods as well. 


Fourier Analysis 

We examine the eigenvalues of the Fourier transform 
of equation (7), and the effects that two precondi- 
tioners have on those eigenvalues. For the spatially 
discrete discontinuous Galerkin method given in equa- 
tion (5) but specialized to one dimension, equation (7) 
becomes 


K(\,\Ji) = a|aQ;+ B i+ iQ; +C,_iQ/-i ] 


U i 


and Q/ = AU/ 4* C /+ iU/+i 4 B ; _ ^ U/ - (8) 


The subscript / is a shorthand for the double sub- 
script /,/=fcl, and 4>i±l is assigned as follows: <Pi+l — 1 
<p t _ i = 0. Eliminating Q gives 


tt(A,U,) = A | [a + B j+ ij [AU; 

+ (c, +i u /+l + B M u,j] 

4 i [AU|_i 4- 

+ (c I+ iU / + B I .iU l . 1 )]}-U, = 0 

= AC,_i (A + B,_i) U,_, 

+ {a [(a + B (+ ±) (a + B ( _i) 

+ C,_iC i+ i] -i}u, 

+ A (a + Bj_j_i^ C j+ xU (+1 
= ^(A)U ; _! + B(A)U, + C(A)U, +1 . (9) 

The Fourier transform is obtained by substituting 
U, = Ue iel where i = and 0 < 9 < n, to give 

tt(A,U) = (U(A)e- is + B(A) + C(A)e'") U (J()) 
= h{ A,0)U 

The preconditioned residual has the form P l /v(A, v). 
The preconditioning operator for block Jacobi is ob- 
tained from equation (9) by retaining only the terms 


involving U / . 

P bj = B(A)Ui 

= {a[(a + B (+ i) (a + B,_i) 

+ c i-i c /+i] - i}u/ 

F% = »(A). 

For block Gauss-Seidel, preconditioner is defined by 
dropping terms involving [//+ 1 ■ 

P gs = &(\)Ui -f *4(A)U/_i 

= {a [(a + Bj + i) (a + B/_i) 

+ c i-i c (+i] - Vi 
+ AC,_ i (a + B ( _i) U;_, 

P 9 , = (B(A)e-’ < ' + A(A)) . 

Analysis of Eigenvalues 

The eigenvalues of 7Z and P -1 7£ are computed at 
50 values of 6 that are uniformly distributed between 
0 and 7 T, and for A ranging from 10" 8 to 10 8 . If p 
denotes the degree of the basis functions, then in one 
dimension there are p 4 1 eigenvalues for each value 
of 0. Figure 1 shows the real component of all eigen- 
values of 1Z for p equals 1-5 and with A — 1.0 (the 
imaginary part is zero). The solid line is proportional 
to p' 1 and closely matches the growth in the magni- 
tude of the eigenvalues. Figures 2(a) and 2(b) show 
the eigenvalues when block Jacobi and block Gauss- 
Seidel preconditioning are used, respectively. Not only 
are the eigenvalues bounded for all p but they are 
grouped into two or three regions. Both precondi- 
tioners place one eigenvalue near zero, bloc k Jacobi 
has one eigenvalue near —2.0, and all of the remaining 
eigenvalues are mapped exactly into —1.0. The varia- 
tion in the eigenvalues within each region indicates the 
degree to which the eigenvalue varies with 0. A clear 
trend observed in both cases is that the eigenvalues 
become more tightly grouped as p increases. Gauss- 
Seidel is not a symmetric operator, and consequently, 
the imaginary component of its eigenvalues is not zero. 
However, the imaginary component is small relative to 
the real component, so all following results will exam- 
ine either the real component or the absolute value ol 
the eigenvalues (denoted as the absolute eigenvalue). 

The above results are for A = 1, but similar trends 
are observed for large and small A. Figures 3(a) and 
3(b) show the eigenvalues for block Jacobi precondi- 
tioning with A = 100 and 0.01, respectively. With 
A = 100 the eigenvalue pattern is nearly identical to 
that of A = 1. Wit h A = 0.01 the eigenvalues are more 
tightly grouped with respect to 0, but the maximum 
and minimum eigenvalue are not as close to 0.0 and 


(13) 

(14) 


( 11 ) 

( 12 ) 
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—2.0. respectively. Similar trends were observed for 
block Gauss- Seidel. 

While the largest absolute eigenvalue limits the time 
step of an explicit method, the smallest absolute eigen- 
value is the least damped and will limit the conver- 
gence of an iterative method. Clearly both precondi- 
tioners reduce the maximum absolute eigenvalue, but 
this is of little value if the minimum absolute eigen- 
value is equally reduced. Thus, another metric for a 
preconditioner is the ratio of the maximum absolute 
eigenvalue to the minimum absolute eigenvalue. Fig- 
ures 4(a-c) show this ratio of absolute eigenvalues with 
no preconditioning, with block Jacobi precondit ioning, 
and with block Gauss-Seidel preconditioning, respec- 
tively. (Note: with no preconditioning the minimum 
absolute eigenvalue is always 1.0, so the ratio is equal 
to the maximum absolute eigenvalue, and the data is 
just the negative of that shown in figure 1; however, 
the data is redrawn on a log scale to facilitate compari- 
son with the other two methods.) Block Jacobi reduces 
the ratio by about an order of magnitude (at p — 5). 
Block Gauss-Seidel further reduces the ratio by about 
half an order (again at p = 5). So both precondition- 
ers provide considerable improvement in convergence, 
but not the three orders of magnitude that, one might 
infer simply by looking at the reduction of maximum 
eigenvalues. 

Optimal Relaxation Schemes 

The tight grouping of eigenvalues means entire 
groups of the error modes can be rapidly damped, 
sometimes completely damped, by the following simple 
update method 

U k+i = U* + wP“ 1 7fc(A,U*). (15) 

The stability of this relaxation step is given by 
G(u,<t{0)) = |1 + u<t(9) | where <r(0) is any eigen- 
value of the operator P -1 7£. By choosing uj = 1, 
all the error modes with an eigenvalue of —1.0 are 
completely damped. Similarly choosing of — 0.5 will 
rapidly damp the error mode with an eigenvalue near 
-2.0. This leaves a single error mode associated with 
the eigenvalue near zero, which is not readily damped 
by any choice of uj that, is not, unstable for the other 
error modes. The block Gauss-Seidel met, hod can be 
overrelaxed with approaching 2 for large p, but, this 
provides only a slight improvement in the damping 
rate of the slow error mode. 

However, if one considers the combined effect of a 
sequence of relaxation steps 

u fc+1 = U k +LJ k P- 1 n(\,U k )< h = 1 ...A (16) 

each with a different, uq then it is possible to choose 
a sequence of that will rapidly damp all eigen- 
modes. In the following sections, we present several 
criteria for selecting an optimal sequence of for block 


Gauss-Seidel and block Jacobi for a fixed A. We then 
construct, a curve fit that accounts for the dependence 
of the uJk sequence on A. 

Block Gauss-Seidel 

For block Gauss-Seidel, we consider a two step se- 
quence whose stability is given by 

G gs {u>2, (t{0)) = G(1.0, c(0))G(u;2,0-(0)), 

and we choose to obtain an optimal damping rate. 
By choosing uq - 1.0, the first step completely damps 
all error modes with an eigenvalue of —1.0 and the 
analysis can focus the single remaining error mode 
whose eigenvalue is near 0.0. This slowly decaying 
error mode, denoted as <7i(0), is dominated by its real 
component which varies monotonically with 6 . This 
behavior is typical for all p and A examined: p < 5 
and I(T 8 < A < 10 s . 

In the following we examine two criteria for choosing 
The first approach minimizes G gs across the entire 
spectrum of 0 by requiring 

G g8 (u2,<r\( 0)) = Ggs{u> 2,n(*r))- (17) 

Figure 5 shows G gs for cr i (0) (solid) and 
(dashed) as a function of uj 2 for P = 3, and for 
A 0.01.0.1,1.0 and 10.0. The symbols show 
G g *(v 2 , <7i(0)) h)r a uniform sampling of 0 which gives 
an indication of the variation of G g9 with 0 at a fixed 
u/ 2 - The symbols also verify that, eit her Gg S (^y ^i(^)) 
or G ff5 (u/ 2 , cr \ ( 7r) ) bounds G gs (u) 2 ,6) above when A 
is not small, and their intersection defines an opti- 
mal u Jn over all 0. For extremely small A, neither 
G^(cJ 2 ,<xi( 0 )) or G fla (u;o,<ri(7r)) bounds G gt ; however 
equation (17) is still a suitable criterion for minimiz- 
ing G gs over all 0. Although the the damping rate is 
stable for all A. the damping rate rises quickly with A, 
and is unacceptablely high at A = 10.0. Also, when 
A is large, G gs (u/ 2 , <Ti(7r)) grows rapidly with and 
quickly becomes unstable. The sensitivity to uq is an 
issue because t he analysis can only be performed for 
the idealized case of a uniform grid. In a realistic case 
in which A varies due to non-uniformity in the mesh, 
the actual optimal value of u 2 may be different from 
that predicted by the analysis, and the predicted op- 
timal value may result in an unstable method. 

A second and more robust approach is to minimize 
Ggs in the high frequency portion of the spectrum (0 > 
k/2) to obtain a relaxation scheme that is well suited 
as a multigrid ' smoother. To this end. we require that 

Ggs(uJ 2,C r l(7r/2)) = cr i( 7r ))- 0®) 

Figure 6(a) shows the u/ 2 obtained by solving equation 
(18) for a range of A. Also shown is the curve defined 
bv 


uJ m + 1 + (u ; m - 1) tanh(« * log io(A) + b) 
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which provides an accurate fit ofu/’o for all A. The coef- 
ficients a and b are obtained from a least square fit over 
the linear subrange seen in figure 6(b), and uj m is the 
maximum value of taken here to be the value of jJo 
when A = 10 s . From figure 7, which shows the depen- 
dence of G gs on u> 2 at, A = 10 6 (similar to figure 5), we 
can see that large variations in u ,' 2 can occur before the 
method becomes unstable, and that small variations 
do not significantly degrade the damping rate. Hence, 
this criterion produces relaxation methods that are ro- 
bust with regard to any discrepancies that might occur 
between the predicted optimal us? and actual values. 
Figures 8(a) and 8(b) show the behavior of G g $ (uj / ) 
as a function of 0 and A, respectively. The maximum 
value over the high frequency spectrum generally oc- 
curs at 9 — 7T /2, and rises to 0.53 for large A. So 
one would expect this relaxation method to be an ex- 
cellent smoother for a multigrid method. It should 
be noted that, because the first, relaxation step com- 
pletely damps all error modes of this idealized linear 
problem except for the cq(0) mode, it is possible to 
perform a single relaxation step w’ith u? = 1.0 followed 
by any number of relaxation steps with lj = . In this 

case the asymptotic convergence rate would reduce to 
(0.53) 2 = 0.28. In the above discussion, p = 3 has 
been used to illustrate the procedure. The same anal- 
ysis has been applied to all p < 5 and the results are 
summarized in Table 1. We find that the asymptotic 
convergence rate is essentially constant as p increases. 

Block Jacobi 

For block Jacobi, we consider a sequence of three 
relaxation steps in which w\ = 0.5, W2 — 1.0 and 
is chosen to provide optimal damping of the high 
frequency portion of the spectrum, 9 > tt/ 2 . The sta- 
bility of such a method is given by 

G bj (u>3 9 <r(9)) = T/(0.5, <r(0))G(l.O, <t(0))G(uj3, (?{&)) 

(20) 

where a(0) is any eigenvalue of P b j [ 7i. Also let ^(0) 
and a 2 (0) denote the maximum and minimum eigen- 
values, respectively. Because the the minimum eigen- 
value <72(0) is near but not, exactly —2.0, the associated 
error mode is not completely damped by the relax- 
ation step with uq = 0.5 and this eigenmode must 
be considered when determining an optimal U3. In 
particular, we choose u/3 by examining the intersec- 
tions of the functions Gbj(^ 3 i ^(tt)), <^ 2 ( ^ ) ) ? 

G bj (u 3 , tri ( tt/ 2)) and G bj {u 3 ,<T'>(n/2)). As shown in 
figure 9 for p — 3 and A = 1.0, the optimal value 
of lu ’3 is associated with the maximum intersection of 
these curves. For the case shown in figure 9, <^3 & 6.7. 
This process is repeated for A ranging from 10 8 to 
10 8 and fitted by equation (19) as done previously for 
block Gauss- Seidel. Figures 10(a) and 10(b) show the 
behavior of 

G m (u>3, 9) = (max(G6j(w3. ^i(0)), Gfcj (u/3, ^(fl)))) 1 ^ 1 


with respect, to changes in 0 and W3, respectively. The 
damping rate for the high-frequency portion of the 
spectrum does not increase above ^ 0.76 as A be- 
comes large. Also, the damping rate increases slowly 
as u.’3 varies from the optimal value indicating that 
the method is robust in this respect. This analysis has 
been performed for all p < 5 and results are summa- 
rized in 'fable 2. As seen previously, the asymptotic 
convergence rate is nearly constant with respect to p . 

Numerical Examples 

In this section we provide a few numerical exam- 
ples that show' that the performance predicted by the 
analysis are realistic and achievable. Details of the 
numerical implementations are not, provided here but 
are throughly described in the references 3 and 7. The 
block Jacobi preconditioned discontinuous Galerkin 
method is applied to the one dimensional heat, equa- 
tion with periodic boundary conditions 

du/dt = V 2 n 

u( 0, .r) = sin(7r,r) - 1 < x < 1 

in which the solution is approximated by polynomi- 
als of degree p = 5. Tests are performed in which 
the solution is evolved in time using either backward 
Euler or third order implicit Runge-Kutta and with 
At /(Ax) 2 = 10. How r ever, here we are interested only 
the convergence of the temporal residual at each time 
step, or each Runge-Kutta stage, which is similar for 
both evolution methods. Figure 1 1 shows the conver- 
gence histories of a typical stage obtained by block 
Jacobi relaxation with constant, lj — 1.0 and opti- 
mal with the three-stage u sequence. Also shown is 
the convergence of block Jacobi relaxation with opti- 
mal three stage uJ sequence combined with a two-level 
multigrid method. The multigrid method is a stan- 
dard V-cycle with 6 relaxation steps on the fine grid 
between coarse grid cycles. The coarse grid solution 
is converged to machine zero by block Jacobi relax- 
ation to simulate the effect, of multiple grid levels. 
The case using the optimal u sequence without multi- 
grid shows some improvement in the clamping rate 
over that of constant u>. The number of relaxation 
steps required to drive the residual to 10“ 12 drops 
from 9706 to 1731; however, this improvement is not 
sufficient to make this approach practical. With multi- 
grid, the convergence improves dramatically and only 
276 relaxation steps are required to drive the resid- 
ual to 10 -12 . This corresponds to a convergence rate 
of IQ " 12 / 276 = 0.905, and represents a 25 fold reduc- 
tion in work when compared to a sixth-order explicit, 
Runge-Kutta method which has an explicit stability 
constraint 3 of At/ A x 2 = 0.0015. 

Extension to Two Dimensions 

The following sections describe the results of the two 
dimensional analysis, construction of optimal relax- 
ation schemes for block Jacobi and numerical results. 
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Two Dimensional Analysis 

In two dimensions, equation (7) has t he form 

H(yUij) = >l(A)U i _ 1 j+B(A)U i j+C(A)U I - +1 j 
4- X>(A)Uij+i 4- (21) 

The Fourier transform is obtained by substituting 

Ui j — Ue*^ r,+ * y ^ to give 

^(A.Uij) = A(\)e-' i9 * +B(\)+C(\)e ie * 

+ T>{\)e ie * + £(A)e _,l9v 
= -h(\,0r,8y)U ( 22 ) 

The preconditioning operator for block Jacobi is sim- 

ply Pbj — B(A)U/ and P bj — #(A). 

Typical eigenvalues are shown in figures 12(a) - 
12(c). There are 21 eigenvalues for the fifth-order 
method (p ~ 4). Figures 12(a) and 12(b) show the 
range of each eigenvalue as the wave numbers vary 
over the complete range 0 < O^.Oy < 7T. Figure 12(c) 
shows a contour plot of the fifth eigenvalue. The anal- 
ysis predicts that both block Gauss-Seidel and block 
Jacobi preconditioners bound the eigenvalues, as seen 
in the one dimensional analysis. Several eigenvalues 
are mapped into - 1 . 0 , however, the remaining eigen- 
values are not tightly grouped but are distributed over 
the range 0 to —1.0 for block Gauss-Seidel and 0 to 
— 2.0 for block Jacobi. 

O ptimal Relaxation 

The eigenvalues are not tightly grouped making it 
necessary to design a relaxation scheme that optimally 
damps a broad spectrum. We choose a sequence of 
three relaxation steps, similar to the previous approach 
for block Jacobi; however, all three relaxation factors 
are allowed to vary freely. The amplification is given 
by 

i ^2 , rr) = Cj (uq , <t)Cj (u/’o , & )G (W 3 , <x) (23) 

The objective of the optimization, illustrat ed in figure 
13, is to minimize G->d f° r ah <r in the range — 2.0 < 
(T < c * . <J * denotes the largest eigenvalue in the high- 
frequency portion of the spectrum and it varies with 
A and p. Note that we are not explicitly concerned 
with the dependence of a on 0. This optimization is 
accomplished by equating 

G 2c/({T* ) = Gndi&a) — ( 6 7 2d(&b) — G‘jrf(~2. 0) (24) 

where <r a and denote a at the local extrema of 
G f 2 d(cr) and are given by dGod/da- = 0.0. Equation 
(24) is quadratic, which allows cr a and <r 6 to be de- 
termined analytically as functions of uq, wo and LO 3 . 
Equation (24) can now be solved to determine uq , u,'o 
and icq . Table 3 gives a * and the optimal w sequence 
for a range of A and for p = 4. uq and uq> are nearly 


constant except at extremely small A. and u / 3 is accu- 
rately fitted by a tanh distribution to give 

u ; 3 = 3.72 4 - 2.72 tanh ( 1 .022 log 10 (A) 4 1 .469) 

Figure 14 shows the asymptotic convergence rate pre- 
dicted by the analysis which approaches 0.93 as A 
becomes large. 

Numerical E xamples 

Figure 15 shows the convergence history from sev- 
eral numerical tests. The test case is the two- 
dimensional heat equation with periodic boundary 
conditions and with the same initial condition as used 
in the one dimensional test described previously. All 
results are for the fifth-order (p = 4) discontinuous 
Galerkin method with block Jacobi preconditioning 
and with At/ Ar = 10. As seen in one dimension, use 
of the optimal u! sequence without multigrid provides 
only a small improvement over the constant uj = 1.0 
case. When relaxation with the optimal u sequence is 
combined with multigrid, the solution converges at a 
spectral radius of 0.89. This rate is faster than that 
predicted by the analysis and represents a 14 fold re- 
duction in work as compared to an explicit method. 

Con el uding Remarks 

Block Jacobi and block Gauss-Seidel relaxation 
methods have been analyzed for the solution ot the 
discontinuous Galerkin method applied to diffusion. 
A Fourier eigenvalue analysis indicates that both pre- 
conditioners bound the eigenvalues. In one dimen- 
sion. the eigenvalues are grouped into well defined 
regions, both preconditioners place one eigenvalue 
near zero, block Jacobi has one eigenvalue near -2, 
and all of the remaining eigenvalues are mapped ex- 
actly into —1.0. In two dimension, several eigenvalues 
are mapped into — 1 . 0 ; however, the remaining eigen- 
values are distributed within the bounding region. In 
general, the eigenvalues become more tightly grouped 
and the minimum eigenvalue moves closer to zero as p 
increases. Block Jacobi reduces the ratio of the maxi- 
mum eigenvalue to the minimum eigenvalue by about 
an order of magnitude; block Gauss-Seidel reduces the 
ratio by an additional half an order. Optimal relax- 
ation methods that employ a sequence of under- and 
over-relaxation are formulated that provided bounded 
convergence rates of the high frequencies as the time 
step increases towards infinity. Numerical examples 
with block Jacobi and multigrid demonstrated conver- 
gence rates near the predicted values. When the use of 
large time steps is appropriate, or when steady state 
solutions are desired, the methods described here of- 
fer more than an order of magnitude reduction in the 
work. 

References 

1 Harold Atkins, Abdelkader Baggag. Can Ozturan, and 
David Keyes, “Parallelization of an Object-Oriented Unstruc- 


6 

American Institute of Aeronautics anl> Astronautics 


AIAA-01-2554 


turcd Aeroacoustic Solver,” Ninth SIAM Conference on Parallel 
Processing for Scientific Computing, March 22-24. 1999. 

2 Harold Atkins, Chi- Wang Shu, “Quadrature-Free Imple- 
mentation of Discontinuous Galerkin Method for Hyperbolic 
Equations” , A/AA Journal Vol. 36, No. 5, 1998, pp. 775-778. 

3 Harold Atkins, Chi- Wang Shu, “Analysis of the Discon- 
tinuous Galerkin Method Applied to the Diffusion Operator 1 , 
AIAA Paper 99-3306, 1999. 

4 J. Tinsley Oden, Ivo Babuska, and Carlos Erik Baumann, 
“A Discontinuous hp Finite Element Method for Diffusion Prob- 
lems,” Journal of Computational Physics , Vol. 146, 1998, pp. 
1-29. 

5 Bernardo Cockburn, Chi- Wang Shu, “T he Local Discon- 
tinuous Galerkin Method For Time-Dependent Convection- 
Diffusion Systems,” SIAM Journal of Numerical Analysis , \ol. 
35, 1998, pp. 2440-2463. 

6 Harold Atkins, David Lockard “A High-Order method us- 
ing Unstructured Grids for the Aeroacoustic Analysis of Realis- 
tic Aircraft Configurations”, AIAA Paper 99-1945, 1999. 

7 Brandt, A., “Multi-level adaptive solutions to boundary- 
value problems.” Mathematics of Computation, Vol 31. pp. 333- 
390, 1977. 



Fig. 1 Eigenvalues of residual operator, A = 1.0 


Table 1 u >/ curve fit parameters and asymptotic 
convergence rate for block Gauss-Seidel 


p 


a 

b 

sjG gs (u> h a{Tc/2)) 
for A = I0 8 

1 

2.32 

1.1858 

0.314 

0 427 

2 

3.85 

1.1486 

0.623 

0.551 

3 

6.15 

1.1144 

0.777 

0.534 

4 

9.13 

1.0878 

0.919 

0.549 

5 

12.79 

1.0781 

0.917 : 

0.557 


Table 2 w/ curve fit parameters and asymptotic 
convergence rate for block Jacobi 


P 


a 

b 

y/G a „(uif,(r(n/2)) 
for A = 10 8 

1 

2.77 

1 . 1919 

0.218 

0.684 

2 

5.26 

1.1379 

0.451 

0.740 

3 

8.76 

1.1 194 

0.556 

0.762 

4 

13.25 

1.0699 

0.495 

0.773 

5 

18.75 

1.0200 

0.602 

0.779 


Table 3 <r* and optimal w for block Jacobi precon- 

ditioning with p = 4 


A 


u/’i 

id r 2 

U/;j 

10~ 2 

-0.297 

0.5302 

0.87057 

2.4309 

0.1 

-0.0845 

0.5343 

0.95946 

4.6989 

1.0 

-0.0308 

0.5353 

0.98483 

6.1458 

10.0 

-0.0237 

0.5354 

0.98830 

6.4078 

10 8 

-0.0228 

0.5355 

0.98873 

6.4413 


(a) 


Re(o) 



(b) 


Re(a) 



Fig. 2 Eigenvalues of preconditioned residual op- 
erator with A = 1.0 and using a) block Jacobi 
preconditioning and b) block Gauss-Seidel precon- 
ditioning. 
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Fig. 3 Eigenvalues of preconditioned residual op- 
erator with block Jacobi preconditioning and with: 
a) A = 100.0 and b) A = 0.01. 



Fig. 4 Ratio of maximum and minimum eigenval- 
ues with A = 1.0 and a) no preconditioning b) block 
Jacobi preconditioning and c) block Gauss-Seidel 
preconditioning. 



Fig. 5 Optimal w 2 over all 0 located at intersection 
of G ffJ (a;2,flr(0)) and G gs (u 2 , <7(tt)). 



(b) 


4 

2 

-2 

-4!- 


OOO O 


9 0 9 


Fig. 6 a) Optimal w for all 6 and curve fit, b) 
linear region for curve fit seen in tanh' " ! (2 * u)f — 

U> T 7i X l ))• 



Fig. 7 Optimal u >2 for all 0 > it/'l located at inter- 
section of G gs (jj 2 , <r ( tt / 2 ) ) and G g ,<{^' 2 , <r( tt) ). 
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Fig. 8 Amplification of relaxation with block 
Gauss-Seidel preconditioning using optimal us se- 
quence: a) (GY*)* VS 0 for A = 0 .1, LO, 10 and 100, 
b) (GY ) * at 9 > 7 r/2. 



5 10 15 


Fig. 10 Amplification of relaxation with block 
Jacobi preconditioning using optimal ,*> sequence, 
Cm = (G b j) 1/3 : a) G m Vs 8 for A =0.1. 1.0, 10 and 
100, b) G f m insensitive to variation in 



Fig. 9 Optimal located at the maximum inter- 
section. 



iteration 

Fig. 11 Convergence of block Jacobi for heat equa- 
tion in one dimension: a) constant u> = 1.0 without 
multigrid, b) optimal jj sequence without multi- 
grid, c) optimal u ; sequence with multigrid. 
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eigenmode, k 




Fig. 14 Predicted asymptotic convergence rate for 
block Jacobi in two dimensions using the optimal 
lj sequence with multigrid, p = 4. 


Fig. 12 Eigenvalues of block Gauss-Seidel and 
block Jacobi in two dimensions: a) distribution 

of all eigenmodes for block Gauss-Seidel precondi- 
tioning, b) distribution of all eigenmodes for block 
Jacobi preconditioning, c) Fourier distribution of 
k = 5 eigenmode for block Gauss-Seidel. 



Fig. 13 Optimal amplification for all rr such that 

<T* < (T < 2 . 0 . 



iteration 

Fig. 15 Convergence of block Jacobi for heat equa- 
tion in two dimensions: a) constant u — 1.0 without 
multigrid, b) optimal w sequence without multi- 
grid, c) constant jj = 1.0 with multigrid, d) optimal 
u> sequence with multigrid. 
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